/*
 *  (c) 2004 Iowa State University
 *      see the LICENSE file in the top level directory
 */

/*	RysPolynomial.c
	Adapted version of the Rys polynomial code in GAMESS.

	converted from GAMESS FORTRAN source 9-1997 BMB

	Routines:

*/
#include <math.h>
#include "Globals.h"
#include "RysPolynomial.h"

#define r12		.275255128608411
#define pie4	.785398163397448
#define r22		2.72474487139158
#define w22		.0917517095361369
#define r13		.190163509193487
#define r23		1.78449274854325
#define w23		.177231492083829
#define r33		5.52534374226326
#define w33		.00511156880411248


void Root123(Root * ldata) {

    double e, y, a1, a2, f1, f2, t1, t2, t3;

/*             *****   VERSION FEBRUARY 13,1975   ***** */
	double x = ldata->x;
	if (x <= 5.0) {
		if (x <= 3.e-7) {	// X IS APPROXIMATELY ZERO.
			switch (ldata->nRoots) {
				case 1:
					ldata->root[0] = 0.5 - x/5.0;
					ldata->weight[0] = 1.0 - x/3.0;
				break;
				case 2:
					ldata->root[0] = .130693606237085 - x * .0290430236082028;
					ldata->root[1] = 2.86930639376291 - x * .637623643058102;
					ldata->weight[0] = .652145154862545 - x * .122713621927067;
					ldata->weight[1] = .347854845137453 - x * .210619711404725;
				break;
				case 3:
					ldata->root[0] = .0603769246832797 - x * .00928875764357368;
					ldata->root[1] = .776823355931043 - x * .119511285527878;
					ldata->root[2] = 6.66279971938567 - x * 1.02504611068957;
					ldata->weight[0] = .467913934572691 - x * .0564876917232519;
					ldata->weight[1] = .360761573048137 - x * .149077186455208;
					ldata->weight[2] = .171324492379169 - x * .127768455150979;
				break;
			}
		} else if (x <= 1.0) {	// x is between 0 and 1
			if (ldata->nRoots <= 2) {
				f1 = ((((((((x * -8.36313918003957e-8 + 1.21222603512827e-6) * 
					x - 1.15662609053481e-5) * x + 9.25197374512647e-5) *
					x - 6.40994113129432e-4) * x + .00378787044215009) *
					x - .0185185172458485) * x + .0714285713298222) *
					x - .199999999997023) * x + .333333333333318;
				ldata->weight[0] = (x + x) * f1 + exp(-x);
			}
			switch (ldata->nRoots) {
				case 1:
					ldata->root[0] = f1 / (ldata->weight[0] - f1);
				break;
				case 2:
					ldata->root[0] = (((((((x * -2.35234358048491e-9 + 2.49173650389842e-8) * 
						x - 4.558315364581e-8) * x - 2.447252174587e-6) * 
						x + 4.743292959463e-5) * x - 5.33184749432408e-4) * 
						x + .00444654947116579) * x - .0290430236084697) * 
						x + .130693606237085;
					ldata->root[1] = (((((((x * -2.4740490232917e-8 + 2.36809910635906e-7) * 
						x + 1.83536773631e-6) * x - 2.066168802076e-5) * 
						x - 1.345693393936e-4) * x - 5.88154362858038e-5) * 
						x + .0532735082098139) * x - .637623643056745) * 
						x + 2.86930639376289;
					ldata->weight[1] = ((f1 - ldata->weight[0]) * ldata->root[0] + f1) *
						(1.0 + ldata->root[1])/(ldata->root[1] - ldata->root[0]);
					ldata->weight[0] -= ldata->weight[1];
				break;
				case 3:
					ldata->root[0] = ((((((x * -5.1018669153887e-10 + 2.4013441570345e-8) * 
						x - 5.01081057744427e-7) * x + 7.58291285499256e-6) *
						x - 9.55085533670919e-5) * x + .00102893039315878) *
						x - .00928875764374337) * x + .060376924683281;
					ldata->root[1] = ((((((x * -1.29646524960555e-8 + 7.74602292865683e-8) * 
						x + 1.56022811158727e-6) * x - 1.58051990661661e-5) *
						x - 3.30447806384059e-4) * x + .00974266885190267) *
						x - .119511285526388) * x + .776823355931033;
					ldata->root[2] = ((((((x * -9.28536484109606e-9 - 3.02786290067014e-7) * 
						x - 2.507344770642e-6) * x - 7.32728109752881e-6) * 
						x + 2.44217481700129e-4) * x + .0494758452357327) * 
						x - 1.02504611065774) * x + 6.66279971938553;
					f2 = ((((((((x * -7.6091148609885e-8 + 1.09552870123182e-6) * 
						x - 1.03463270693454e-5) * x + 8.16324851790106e-5) *
						x - 5.55526624875562e-4) * x + .00320512054753924) *
						x - .015151513983854) * x + .0555555554649585) * 
						x - .142857142854412) * x + .199999999999986;
					e = exp(-x);
					f1 = ((x + x) * f2 + e) / 3.0;
					ldata->weight[0] = (x + x) * f1 + e;
					t1 = ldata->root[0] / (ldata->root[0] + 1.0);
					t2 = ldata->root[1] / (ldata->root[1] + 1.0);
					t3 = ldata->root[2] / (ldata->root[2] + 1.0);
					a2 = f2 - t1 * f1;
					a1 = f1 - t1 * ldata->weight[0];
					ldata->weight[2] = (a2 - t2 * a1) / ((t3 - t2) * (t3 - t1));
					ldata->weight[1] = (t3 * a1 - a2) / ((t3 - t2) * (t2 - t1));
					ldata->weight[0] -= (ldata->weight[1] + ldata->weight[2]);
				break;
			}
		} else if (x <= 3.0) {	//x is between 1 and 3
			y = x - 2.0;
			if (ldata->nRoots <= 2) {
				f1 = ((((((((((y * -1.61702782425558e-10 + 1.96215250865776e-9) * y -
					2.14234468198419e-8) * y + 2.17216556336318e-7) * y -
					1.98850171329371e-6) * y + 1.62429321438911e-5) * y -
					1.16740298039895e-4) * y + 7.24888732052332e-4) * y -
					.00379490003707156) * y + .0161723488664661) * y -
					.0529428148329736) * y + .115702180856167;
				ldata->weight[0] = (x + x) * f1 + exp(-x);
			}
			switch (ldata->nRoots) {
				case 1:
					ldata->root[0] = f1 / (ldata->weight[0] - f1);
				break;
				case 2:
					ldata->root[0] = (((((((((y * -6.36859636616415e-12 + 8.4741706477627e-11) * y -
						5.152207846962e-10) * y - 3.846389873308e-10) * y +
						8.47225338838e-8) * y - 1.85306035634293e-6) * y +
						2.47191693238413e-5) * y - 2.49018321709815e-4) * y +
						.00219173220020161) * y - .0163329339286794) * y + .0868085688285261;
					ldata->root[1] = (((((((((y * 1.45331350488343e-10 + 2.07111465297976e-9) * y -
						1.878920917404e-8) * y - 1.725838516261e-7) * y +
						2.247389642339e-6) * y + 9.76783813082564e-6) * y -
						1.93160765581969e-4) * y - .00158064140671893) * y +
						.0485928174507904) * y - .430761584997596) * y + 1.8040097453795;
					ldata->weight[1] = ((f1 - ldata->weight[0]) * ldata->root[0] + f1) *
						(1.0 + ldata->root[1])/(ldata->root[1] - ldata->root[0]);
					ldata->weight[0] -= ldata->weight[1];
				break;
				case 3:
					ldata->root[0] = ((((((((y * 1.44687969563318e-12 + 4.85300143926755e-12) * y - 
						6.55098264095516e-10) * y + 1.56592951656828e-8) * y - 
						2.60122498274734e-7) * y + 3.86118485517386e-6) * y - 
						5.13430986707889e-5) * y + 6.03194524398109e-4) * y - 
						.0061121934982509) * y + .0452578254679079;
					ldata->root[1] = (((((((y * 6.95964248788138e-10 - 5.35281831445517e-9) * y - 
						6.745205954533e-8) * y + 1.502366784525e-6) * y + 
						9.923326947376e-7) * y - 3.89147469249594e-4) * y + 
						.00751549330892401) * y - .08487781203634) * y + .573928229597613;
					ldata->root[2] = ((((((((y * -2.81496588401439e-10 + 3.61058041895031e-9) * y + 
						4.53631789436255e-8) * y - 1.40971837780847e-7) * y - 
						6.05865557561067e-6) * y - 5.15964042227127e-5) * y + 
						3.34761560498171e-5) * y + .0504871005319119) * y - 
						.824708946991557) * y + 4.81234667357205;
					f2 = ((((((((((y * -1.4804423107214e-10 + 1.78157031325097e-9) * y - 
						1.92514145088973e-8) * y + 1.92804632038796e-7) * y - 
						1.73806555021045e-6) * y + 1.39195169625425e-5) * y - 
						9.74574633246452e-5) * y + 5.83701488646511e-4) * y - 
						.00289955494844975) * y + .011384700111381) * y - 
						.0323446977320647) * y + .0529428148329709;
					e = exp(-x);
					f1 = ((x + x) * f2 + e) / 3.0;
					ldata->weight[0] = (x + x) * f1 + e;
					t1 = ldata->root[0] / (ldata->root[0] + 1.0);
					t2 = ldata->root[1] / (ldata->root[1] + 1.0);
					t3 = ldata->root[2] / (ldata->root[2] + 1.0);
					a2 = f2 - t1 * f1;
					a1 = f1 - t1 * ldata->weight[0];
					ldata->weight[2] = (a2 - t2 * a1) / ((t3 - t2) * (t3 - t1));
					ldata->weight[1] = (t3 * a1 - a2) / ((t3 - t2) * (t2 - t1));
					ldata->weight[0] -= (ldata->weight[1] + ldata->weight[2]);
				break;
			}
		} else {	//x is between 3 and 5
			y = x - 4.0;
			if (ldata->nRoots <= 2) {
				f1 = ((((((((((y * -2.62453564772299e-11 + 3.24031041623823e-10) * y - 
					3.614965656163e-9) * y + 3.760256799971e-8) * y - 
					3.553558319675e-7) * y + 3.022556449731e-6) * y - 
					2.290098979647e-5) * y + 1.526537461148e-4) * y - 
					8.81947375894379e-4) * y + .00433207949514611) * y - 
					.0175257821619926) * y + .0528406320615584;
				ldata->weight[0] = (x + x) * f1 + exp(-x);
			}
			switch (ldata->nRoots) {
				case 1:
					ldata->root[0] = f1 / (ldata->weight[0] - f1);
				break;
				case 2:
					ldata->root[0] = ((((((((y * -4.11560117487296e-12 + 7.10910223886747e-11) * y - 
						1.73508862390291e-9) * y + 5.93066856324744e-8) * y - 
						9.76085576741771e-7) * y + 1.08484384385679e-5) * y - 
						1.12608004981982e-4) * y + .00116210907653515) * y - 
						.00989572595720351) * y + .0612589701086408;
					ldata->root[1] = (((((((((y * -1.80555625241001e-10 + 5.44072475994123e-10) * y + 
						1.60349804524e-8) * y - 1.497986283037e-7) * y - 
						7.017002532106e-7) * y + 1.85882653064034e-5) * y - 
						2.04685420150802e-5) * y - .00249327728643089) * y + 
						.0356550690684281) * y - .260417417692375) * y + 1.12155283108289;
					ldata->weight[1] = ((f1 - ldata->weight[0]) * ldata->root[0] + f1) *
						(ldata->root[1] + 1.) / (ldata->root[1] - ldata->root[0]);
					ldata->weight[0] -= ldata->weight[1];
				break;
				case 3:
					ldata->root[0] = (((((((y * 1.44265709189601e-11 - 4.66622033006074e-10) * y + 
						7.649155832025e-9) * y - 1.229940017368e-7) * y + 
						2.026002142457e-6) * y - 2.87048671521677e-5) * y + 
						3.70326938096287e-4) * y - .00421006346373634) * y + .0350898470729044;
					ldata->root[1] = ((((((((y * -2.65526039155651e-11 + 1.97549041402552e-10) * y + 
						2.15971131403034e-9) * y - 7.95045680685193e-8) * y + 
						5.15021914287057e-7) * y + 1.11788717230514e-5) * y - 
						3.33739312603632e-4) * y + .00530601428208358) * y - 
						.0593483267268959) * y + .431180523260239;
					ldata->root[2] = ((((((((y * -3.92833750584041e-10 - 4.1642322978228e-9) * y + 
						4.42413039572867e-8) * y + 6.40574545989551e-7) * y - 
						3.05512456576552e-6) * y - 1.05296443527943e-4) * y - 
						6.14120969315617e-4) * y + .0489665802767005) * y - 
						.624498381002855) * y + 3.36412312243724;
					f2 = ((((((((((y * -2.36788772599074e-11 + 2.89147476459092e-10) * y - 
						3.18111322308846e-9) * y + 3.25336816562485e-8) * y - 
						3.00873821471489e-7) * y + 2.48749160874431e-6) * y - 
						1.81353179793672e-5) * y + 1.14504948737066e-4) * y - 
						6.10614987696677e-4) * y + .00264584212770942) * y - 
						.00866415899015349) * y + .0175257821619922;
					e = exp(-x);
					f1 = ((x + x) * f2 + e) / 3.0;
					ldata->weight[0] = (x + x) * f1 + e;
					t1 = ldata->root[0] / (ldata->root[0] + 1.0);
					t2 = ldata->root[1] / (ldata->root[1] + 1.0);
					t3 = ldata->root[2] / (ldata->root[2] + 1.0);
					a2 = f2 - t1 * f1;
					a1 = f1 - t1 * ldata->weight[0];
					ldata->weight[2] = (a2 - t2 * a1) / ((t3 - t2) * (t3 - t1));
					ldata->weight[1] = (t3 * a1 - a2) / ((t3 - t2) * (t2 - t1));
					ldata->weight[0] -= (ldata->weight[1] + ldata->weight[2]);
				break;
			}
		}
	} else {
		if (x <= 10.0) {
    		e = exp(-x);
			ldata->weight[0] = ((((((.46897511375022 / x - .69955602298985) / x + 
				.53689283271887) / x - .32883030418398) / x + .24645596956002) / x -
				.49984072848436) / x - 3.1501078774085e-6) * e + sqrt(pie4 / x);
			f1 = (ldata->weight[0] - e) / (x + x);
			switch (ldata->nRoots) {
				case 1:
					ldata->root[0] = f1 / (ldata->weight[0] - f1);
				break;
				case 2:
					y = x - 7.5;
					ldata->root[0] = (((((((((((((y * -1.43632730148572e-16 + 2.38198922570405e-16) * y +
						1.3583196188e-14) * y - 7.064522786879e-14) * y - 7.719300212748e-13) * y +
						7.802544789997e-12) * y + 6.628721099436e-11) * y - 1.775564159743e-9) * y + 
						1.71382882399e-8) * y - 1.497500187053e-7) * y + 2.283485114279e-6) * y -
						3.76953869614706e-5) * y + 4.74791204651451e-4) * y - .00460448960876139) * y + 
						.0372458587837249;
					ldata->root[1] = ((((((((((((y * 2.487916227989e-14 - 1.36113510175724e-13) * y - 
						2.224334349799e-12) * y + 4.190559455515e-11) * y - 2.222722579924e-10) * y -
						2.624183464275e-9) * y + 6.128153450169e-8) * y - 4.383376014528e-7) * y - 
						2.4995220023291e-6) * y + 1.0323664788832e-4) * y - .00144614664924989) * y +
						.0135094294917224) * y - .0953478510453887) * y + .54476524568679;
					ldata->weight[1] = ((f1 - ldata->weight[0]) * ldata->root[0] + f1) *
						(ldata->root[1] + 1.) / (ldata->root[1] - ldata->root[0]);
					ldata->weight[0] -= ldata->weight[1];
				break;
				case 3:
					f2 = (f1 + f1 + f1 - e) / (x + x);
					y = x - 7.5;
					ldata->root[0] = (((((((((((y * 5.74429401360115e-16 + 7.11884203790984e-16) * y - 
						6.736701449826e-14) * y - 6.264613873998e-13) * y + 1.31541892704e-11) *
						y - 4.23879635610964e-11) * y + 1.39032379769474e-9) * y -
						4.65449552856856e-8) * y + 7.34609900170759e-7) * y -
						1.08656008854077e-5) * y + 1.77930381549953e-4) * y -
						.00239864911618015) * y + .0239112249488821;
					ldata->root[1] = (((((((((((y * 1.1346409620912e-14 + 6.99375313934242e-15) * y - 
						8.595618132088e-13) * y - 5.293620408757e-12) * y - 2.492175211635e-11) *
						y + 2.73681574882729e-9) * y - 1.06656985608482e-8) * y -
						4.40252529648056e-7) * y + 9.68100917793911e-6) * y -
						1.68211091755327e-4) * y + .00269443611274173) * y - .0323845035189063)
						* y + .275969447451882;
					ldata->root[2] = ((((((((((((y * 6.66339416996191e-15 + 1.84955640200794e-13) * y - 
						1.985141104444e-12) * y - 2.309293727603e-11) * y + 3.917984522103e-10) *
						y + 1.663165279876e-9) * y - 6.205591993923e-8) * y + 8.769581622041e-9) *
						y + 8.97224398620038e-6) * y - 3.14232666170796e-5) * y - .00183917335649633) *
						y + .0351246831672571) * y - .32233505127086) * y + 1.7358283175543;
					f1 = ((x + x) * f2 + e) / 3.0;
					ldata->weight[0] = (x + x) * f1 + e;
					t1 = ldata->root[0] / (ldata->root[0] + 1.0);
					t2 = ldata->root[1] / (ldata->root[1] + 1.0);
					t3 = ldata->root[2] / (ldata->root[2] + 1.0);
					a2 = f2 - t1 * f1;
					a1 = f1 - t1 * ldata->weight[0];
					ldata->weight[2] = (a2 - t2 * a1) / ((t3 - t2) * (t3 - t1));
					ldata->weight[1] = (t3 * a1 - a2) / ((t3 - t2) * (t2 - t1));
					ldata->weight[0] -= (ldata->weight[1] + ldata->weight[2]);
				break;
			}
		} else if (x <= 15.0) {
    		e = exp(-x);
			ldata->weight[0] = (((-.18784686463512 / x + .22991849164985) / x - .49893752514047) /
			x - 2.1916512131607e-5) * e + sqrt(pie4 / x);
			f1 = (ldata->weight[0] - e) / (x + x);
			switch (ldata->nRoots) {
				case 1:
					ldata->root[0] = f1 / (ldata->weight[0] - f1);
				break;
				case 2:
					ldata->root[0] = ((((x * -1.01041157064226e-5 + .00119483054115173) * 
						x - .0673760231824074) * x + 1.25705571069895) * x +
						(((-8576.09422987199 / x + 5910.05939591842) / x - 1708.07677109425) /
						x + 264.536689959503) / x - 23.8570496490846) * e + r12 / (x - r12);
					ldata->root[1] = (((x * 3.39024225137123e-4 - .0934976436343509) * x -
						4.2221648330632) * x + (((-2084.57050986847 / x - 1049.99071905664) /
						x + 339.891508992661) / x - 156.184800325063) / x + 8.00839033297501) *
						e + r22 / (x - r22);
					ldata->weight[1] = ((f1 - ldata->weight[0]) * ldata->root[0] + f1) *
						(ldata->root[1] + 1.) / (ldata->root[1] - ldata->root[0]);
					ldata->weight[0] -= ldata->weight[1];
				break;
				case 3:
					f2 = (f1 + f1 + f1 - e) / (x + x);
					y = x - 12.5;
					ldata->root[0] = (((((((((((y * 4.4213300128309e-16 - 2.77189767070441e-15) * y - 
						4.084026087887e-14) * y + 5.379885121517e-13) * y + 1.882093066702e-12) *
						y - 8.67286219861085e-11) * y + 7.11372337079797e-10) * y -
						3.55578027040563e-9) * y + 1.29454702851936e-7) * y -
						4.14222202791434e-6) * y + 8.04427643593792e-5) * y -
						.00118587782909876) * y + .0153435577063174;
					ldata->root[1] = (((((((((((y * 6.85146742119357e-15 - 1.08257654410279e-14) * y - 
						8.579165965128e-13) * y + 6.642452485783e-12) * y + 4.798806828724e-11) *
						y - 1.13413908163831e-9) * y + 7.08558457182751e-9) * y -
						5.59678576054633e-8) * y + 2.51020389884249e-6) * y -
						6.63678914608681e-5) * y + .00111888323089714) * y - .0145361636398178) *
						y + .165077877454402;
					ldata->root[2] = ((((((((((((y * 3.20622388697743e-15 - 2.73458804864628e-14) * y - 
						3.157134329361e-13) * y + 8.654129268056e-12) * y - 5.625235879301e-11) *
						y - 7.718080513708e-10) * y + 2.064664199164e-8) * y - 1.567725007761e-7) *
						y - 1.57938204115055e-6) * y + 6.27436306915967e-5) * y -
						.00101308723606946) * y + .0113901881430697) * y - .10144965289945) * y + .777203937334739;
					f1 = ((x + x) * f2 + e) / 3.0;
					ldata->weight[0] = (x + x) * f1 + e;
					t1 = ldata->root[0] / (ldata->root[0] + 1.0);
					t2 = ldata->root[1] / (ldata->root[1] + 1.0);
					t3 = ldata->root[2] / (ldata->root[2] + 1.0);
					a2 = f2 - t1 * f1;
					a1 = f1 - t1 * ldata->weight[0];
					ldata->weight[2] = (a2 - t2 * a1) / ((t3 - t2) * (t3 - t1));
					ldata->weight[1] = (t3 * a1 - a2) / ((t3 - t2) * (t2 - t1));
					ldata->weight[0] -= (ldata->weight[1] + ldata->weight[2]);
				break;
			}
		} else if (x <= 33.0) {	// x between 15 to 33
			e = exp(-x);
			ldata->weight[0] = ((.1962326414943 / x - .4969524146449) / x - 6.0156581186481e-5) *
				e + sqrt(pie4 / x);
			f1 = (ldata->weight[0] - e) / (x + x);
			switch (ldata->nRoots) {
				case 1:
					ldata->root[0] = f1 / (ldata->weight[0] - f1);
				break;
				case 2:
					ldata->root[0] = ((((x * -1.14906395546354e-6 + 1.76003409708332e-4) * 
						x - .0171984023644904) * x - .137292644149838) * x + (-47.5742064274859 /
						x + 9.21005186542857) / x - .0231080873898939) * e + r12 / (x - r12);
					ldata->root[1] = (((x * 3.64921633404158e-4 - .0971850973831558) * x 
						- 4.02886174850252) * x + (-135.831002139173 / x - 86.6891724287962) /
						x + 2.98011277766958) * e + r22 / (x - r22);
					ldata->weight[1] = ((f1 - ldata->weight[0]) * ldata->root[0] + f1) *
						(ldata->root[1] + 1.) / (ldata->root[1] - ldata->root[0]);
					ldata->weight[0] -= ldata->weight[1];
				break;
				case 3:
					f2 = (f1 + f1 + f1 - e) / (x + x);
					if (x <= 20.) {
						ldata->root[0] = ((((((x * -2.43270989903742e-6 + 3.57901398988359e-4) * 
							x - .0234112415981143) * x + .781425144913975) * x -
							17.3209218219175) * x + 243.517435690398) * x +
							(-19761.1541576986 / x + 9824.41363463929) / x -
							2079.70687843258) * e + r13 / (x - r13);
						ldata->root[1] = (((((x * -2.62627010965435e-4 + .0349187925428138) * 
							x - 3.0933761873188) * x + 107.037141010778) * x -
							2366.59637247087) * x + ((-2916691.1368102 / x +
							1411295.05262758) / x - 291532.335433779) / x + 33520.2872835409) *
							e + r23 / (x - r23);
						ldata->root[2] = (((((x * 9.31856404738601e-5 - .0287029400759565) * 
							x - .783503697918455) * x - 18.4338896480695) * x + 404.996712650414) *
							x + (-189829.509315154 / x + 51149.8390849158) / x -
							6881.45821789955) * e + r33 / (x - r33);
					} else {
						ldata->root[0] = ((((x * -4.97561537069643e-4 - .0500929599665316) * 
							x + 1.31099142238996) * x - 18.8336409225481) * x -
							660.344754467191 / x + 164.931462413877) * e + r13 / (x - r13);
						ldata->root[1] = ((((x * -.00448218898474906 - .517373211334924) * x +
							11.3691058739678) * x - 165.426392885291) * x - 6309.09125686731 /
							x + 1522.31757709236) * e + r23 / (x - r23);
						ldata->root[2] = ((((x * -.0138368602394293 - 1.77293428863008) * x + 
							17.3639054044562) * x - 357.615122086961) * x - 14573.4701095912 /
							x + 2698.31813951849) * e + r33 / (x - r33);
					}
					f1 = ((x + x) * f2 + e) / 3.0;
					ldata->weight[0] = (x + x) * f1 + e;
					t1 = ldata->root[0] / (ldata->root[0] + 1.0);
					t2 = ldata->root[1] / (ldata->root[1] + 1.0);
					t3 = ldata->root[2] / (ldata->root[2] + 1.0);
					a2 = f2 - t1 * f1;
					a1 = f1 - t1 * ldata->weight[0];
					ldata->weight[2] = (a2 - t2 * a1) / ((t3 - t2) * (t3 - t1));
					ldata->weight[1] = (t3 * a1 - a2) / ((t3 - t2) * (t2 - t1));
					ldata->weight[0] -= (ldata->weight[1] + ldata->weight[2]);
				break;
			}
		} else {	// x > 33
			ldata->weight[0] = sqrt(pie4 / x);
			switch (ldata->nRoots) {
				case 1:
					ldata->root[0] = 0.5 / (x - 0.5);
				break;
				case 2:
					if (x <= 40.0) {
						e = exp(-x);
						ldata->root[0] = (x * -.87894730749888 + 10.9243702330261) * e + r12 / (x - r12);
						ldata->root[1] = (x * -9.28903924275977 + 81.0642367843811) * e + r22 /
							(x - r22);
						ldata->weight[1] = (x * 4.468573893084 - 77.9250653461045) * e + w22 *
							ldata->weight[0];
						ldata->weight[0] -= ldata->weight[1];
					} else {
						ldata->root[0] = r12 / (x - r12);
						ldata->root[1] = r22 / (x - r22);
						ldata->weight[1] = w22 * ldata->weight[0];
						ldata->weight[0] -= ldata->weight[1];
					}
				break;
				case 3:
					if (x <= 47.0) {
						e = exp(-x);
						ldata->root[0] = ((x * -7.39058467995275 + 321.318352526305) * x - 
							3994.33696473658) * e + r13 / (x - r13);
						ldata->root[1] = ((x * -73.8726243906513 + 3135.69966333873) * x - 
							38686.2867311321) * e + r23 / (x - r23);
						ldata->root[2] = ((x * -263.750565461336 + 10441.2168692352) * x - 
							128094.577915394) * e + r33 / (x - r33);
						ldata->weight[2] = (((x * .152258947224714 - 8.30661900042651) * x + 
							192.977367967984) * x - 1677.87926005344) * e + w33 * ldata->weight[0];
						ldata->weight[1] = ((x * 61.5072615497811 - 2919.80647450269) * x + 
							38079.4303087338) * e + w23 * ldata->weight[0];
						ldata->weight[0] -= (ldata->weight[1] + ldata->weight[2]);
					} else {
						ldata->root[0] = r13 / (x - r13);
						ldata->root[1] = r23 / (x - r23);
						ldata->root[2] = r33 / (x - r33);
						ldata->weight[1] = w23 * ldata->weight[0];
						ldata->weight[2] = w33 * ldata->weight[0];
						ldata->weight[0] -= (ldata->weight[1] + ldata->weight[2]);
					}
				break;
			}
		}
	}
}

#define r14		.145303521503316
#define r24		1.33909728812636
#define w24		.234479815323517
#define r34		3.92696350135829
#define w34		.0192704402415764
#define r44		8.58863568901199
#define w44		2.25229076750736e-4

void Root4(Root * ldata) {
	double e, y;

/*          *****   VERSION FEBRUARY 16,1975   ***** */

	double x = ldata->x;
	if (x <= 15.0) {
		if (x <= 3e-7) {	//x is almost zero
			ldata->root[0] = .0348198973061471 - x * .00409645850660395;
			ldata->root[1] = .381567185080042 - x * .0448902570656719;
			ldata->root[2] = 1.73730726945891 - x * .204389090547327;
			ldata->root[3] = 11.8463056481549 - x * 1.39368301742312;
			ldata->weight[0] = .362683783378362 - x * .0313844305713928;
			ldata->weight[1] = .313706645877886 - x * .0898046242557724;
			ldata->weight[2] = .222381034453372 - x * .129314370958973;
			ldata->weight[3] = .101228536290376 - x * .0828299075414321;
		} else if (x <= 1.0) {	// x zero to 1
			ldata->root[0] = ((((((x * -1.95309614628539e-10 + 5.19765728707592e-9) * x -
				1.01756452250573e-7) * x + 1.72365935872131e-6) * x - 2.61203523522184e-5) *
				x + 3.5292130876988e-4) * x - .00409645850658433) * x + .0348198973061469;
			ldata->root[1] = (((((x * -1.89554881382342e-8 + 3.07583114342365e-7) * x +
				1.270981734393e-6) * x - 1.417298563884e-4) * x + .003226979163176) * x -
				.0448902570678178) * x + .381567185080039;
			ldata->root[2] = ((((((x * 1.77280535300416e-9 + 3.36524958870615e-8) * x -
				2.58341529013893e-7) * x - 1.1364489566232e-5) * x - 7.91549618884063e-5) *
				x + .0103825827346828) * x - .204389090525137) * x + 1.73730726945889;
			ldata->root[3] = (((((x * -5.61188882415248e-8 - 2.4948073307246e-7) * x +
				3.428685057114e-6) * x + 1.679007454539e-4) * x + .04722855585715) * x -
				1.39368301737828) * x + 11.8463056481543;
			ldata->weight[0] = ((((((x * -1.14649303201279e-8 + 1.88015570196787e-7) * x -
				2.33305875372323e-6) * x + 2.68880044371597e-5) * x - 2.94268428977387e-4) *
				x + .00306548909776613) * x - .0313844305680096) * x + .362683783378335;
			ldata->weight[1] = ((((((((x * -4.11720483772634e-9 + 6.54963481852134e-8) * x -
				7.20045285129626e-7) * x + 6.93779646721723e-6) * x - 6.05367572016373e-5) * x +
				4.74241566251899e-4) * x - .00326956188125316) * x + .0191883866626681) * x -
				.0898046242565811) * x + .313706645877886;
			ldata->weight[2] = ((((((((x * -3.41688436990215e-8 + 5.07238960340773e-7) * x -
				5.0167562840822e-6) * x + 4.20363420922845e-5) * x - 3.08040221166823e-4) *
				x + .00194431864731239) * x - .0102477820460278) * x + .0428670143840073) * 
				x - .129314370962569) * x + .222381034453369;
			ldata->weight[3] = (((((((((x * 4.99660550769508e-9 - 7.9458596331012e-8) * x +
				8.359072409485e-7) * x - 7.42236921061e-6) * x + 5.76337430816e-5) * x -
				3.86645606718233e-4) * x + .00218417516259781) * x - .00999791027771119) * 
				x + .034879109737737) * x - .0828299075413889) * x + .101228536290376;
		} else if (x <= 5.0) {	// x between 1 and 5
			y = x - 3.;
			ldata->root[0] = (((((((((y * -1.48570633747284e-15 - 1.33273068108777e-13) * y + 
				4.06854369667e-12) * y - 9.163164161821e-11) * y + 2.046819017845e-9) * y -
				4.03076426299031e-8) * y + 7.29407420660149e-7) * y - 1.23118059980833e-5) *
				y + 1.88796581246938e-4) * y - .00253262912046853) * y + .0251198234505021;
			ldata->root[1] = (((((((((y * 1.35830583483312e-13 - 2.29772605964836e-12) * y - 
				3.821500128045e-12) * y + 6.844424214735e-10) * y - 1.048063352259e-8) * y +
				1.50083186233363e-8) * y + 3.48848942324454e-6) * y - 1.08694174399193e-4) *
				y + .00208048885251999) * y - .0291205805373793) * y + .272276489515713;
			ldata->root[2] = (((((((((y * 5.02799392850289e-13 + 1.07461812944084e-11) * y - 
				1.482277886411e-10) * y - 2.153585661215e-9) * y + 3.654087802817e-8) * y +
				5.1592957583012e-7) * y - 9.52388379435709e-6) * y - 2.16552440036426e-4) *
				y + .0090355146956832) * y - .145505469175613) * y + 1.21449092319186;
			ldata->root[3] = (((((((((y * -1.08510370291979e-12 + 6.41492397277798e-11) * y + 
				7.542387436125e-10) * y - 2.213111836647e-9) * y - 1.448228963549e-7) * y -
				1.95670833237101e-6) * y - 1.07481314670844e-5) * y + 1.49335941252765e-4) *
				y + .0487791531990593) * y - 1.10559909038653) * y + 8.0950202861178;
			ldata->weight[0] = ((((((((((y * -4.65801912689961e-14 + 7.586695071068e-13) * y - 
				1.186387548048e-11) * y + 1.862334710665e-10) * y - 2.799399389539e-9) * y +
				4.148972684255e-8) * y - 5.9335680796e-7) * y + 8.168349266115e-6) * y -
				1.08989176177409e-4) * y + .00141357961729531) * y - .0187588361833659) * y + 
				.289898651436026;
			ldata->weight[1] = ((((((((((((y * -1.46345073267549e-14 + 2.25644205432182e-13) * y 
				- 3.116258693847e-12) * y + 4.32190875661e-11) * y - 5.673270062669e-10) * y +
				7.00629596296e-9) * y - 8.120186517e-8) * y + 8.77529464577e-7) * y -
				8.77829235749024e-6) * y + 8.04372147732379e-5) * y - 6.64149238804153e-4) * y +
				.00481181506827225) * y - .0288982669486183) * y + .156247249979288;
			ldata->weight[2] = (((((((((((((y * 9.06812118895365e-15 - 1.40541322766087e-13) * y 
				+ 1.919270015269e-12) * y - 2.60513573901e-11) * y + 3.299685839012e-10) * y -
				3.86354139348735e-9) * y + 4.16265847927498e-8) * y - 4.0946283547147e-7) * y + 
				3.64018881086111e-6) * y - 2.88665153269386e-5) * y + 2.00515819789028e-4) * y -
				.00118791896897934) * y + .00575223633388589) * y - .0209400418772687) * y + 
				.0485368861938873;
			ldata->weight[3] = ((((((((((((((y * -9.74835552342257e-16 + 1.57857099317175e-14) * 
				y - 2.249993780112e-13) * y + 3.173422008953e-12) * y - 4.16115945968e-11) * y +
				5.021343560166e-10) * y - 5.545047534808e-9) * y + 5.554146993491e-8) * y - 
				4.99048696190133e-7) * y + 3.96650392371311e-6) * y - 2.73816413291214e-5) * y +
				1.60106988333186e-4) * y - 7.64560567879592e-4) * y + .00281330044426892) * y - 
				.00716227030134947) * y + .00966077262223353;
		} else if (x <= 10.0) {	//x between 5 and 10
			y = x - 7.5;
			ldata->root[0] = (((((((((y * 4.64217329776215e-15 - 6.27892383644164e-15) * y + 
				3.462236347446e-13) * y - 2.92722935535e-11) * y + 5.090355371676e-10) * y -
				9.97272656345253e-9) * y + 2.37835295639281e-7) * y - 4.60301761310921e-6) *
				y + 8.42824204233222e-5) * y - .00137983082233081) * y + .0166630865869375;
			ldata->root[1] = (((((((((y * 2.93981127919047e-14 + 8.47635639065744e-13) * y - 
				1.446314544774e-11) * y - 6.149155555753e-12) * y + 8.484275604612e-10) * y -
				6.10898827887652e-8) * y + 2.39156093611106e-6) * y - 5.35837089462592e-5) * y + 
				.00100967602595557) * y - .0157769317127372) * y + .174853819464285;
			ldata->root[2] = ((((((((((y * 2.93523563363e-14 - 6.4004177666702e-14) * y - 
				2.695740446312e-12) * y + 1.027082960169e-10) * y - 5.82203865678e-10) * y -
				3.159991002539e-8) * y + 4.327249251331e-7) * y + 4.856768455119e-6) * y - 
				2.54617989427762e-4) * y + .00554843378106589) * y - .0795013029486684) *
				y + .720206142703162;
			ldata->root[3] = (((((((((((y * -1.62212382394553e-14 + 7.68943641360593e-13) * y + 
				5.764015756615e-12) * y - 1.380635298784e-10) * y - 1.476849808675e-9) * y +
				1.84347052385605e-8) * y + 3.34382940759405e-7) * y - 1.39428366421645e-6) *
				y - 7.50249313713996e-5) * y - 6.26495899187507e-4) * y + .0469716410901162) *
				y - .666871297428209) * y + 4.11207530217806;
			ldata->weight[0] = ((((((((((y * -1.65995045235997e-15 + 6.91838935879598e-14) * y - 
				9.131223418888e-13) * y + 1.403341829454e-11) * y - 3.672235069444e-10) * y +
				6.36696254699e-9) * y - 1.039220021671e-7) * y + 1.959098751715e-6) * y - 
				3.33474893152939e-5) * y + 5.72164211151013e-4) * y - .0105583210553392) *
				y + .226696066029591;
			ldata->weight[1] = ((((((((((((y * -3.57248951192047e-16 + 6.25708409149331e-15) * y 
				- 9.657033089714e-14) * y + 1.507864898748e-12) * y - 2.33252225611e-11) * y +
				3.428545616603e-10) * y - 4.698730937661e-9) * y + 6.21997763513e-8) * y - 
				7.83008889613661e-7) * y + 9.08621687041567e-6) * y - 9.86368311253873e-5) * y
				+ 9.69632496710088e-4) * y - .00814594214284187) * y + .0850218447733457;
			ldata->weight[2] = (((((((((((((y * 1.64742458534277e-16 - 2.6851226592841e-15) * y + 
				3.788890667676e-14) * y - 5.508918529823e-13) * y + 7.555896810069e-12) * y -
				9.69039768312637e-11) * y + 1.16034263529672e-9) * y - 1.28771698573873e-8) *
				y + 1.31949431805798e-7) * y - 1.23673915616005e-6) * y + 1.04189803544936e-5) *
				y - 7.79566003744742e-5) * y + 5.03162624754434e-4) * y - .00255138844587555) *
				y + .0113250730954014;
			ldata->weight[3] = ((((((((((((((y * -1.55714130075679e-17 + 2.57193722698891e-16) * 
				y - 3.626606654097e-15) * y + 5.234734676175e-14) * y - 7.067105402134e-13) *
				y + 8.79351266489e-12) * y - 1.006088923498e-10) * y + 1.050565098393e-9) * y - 
				9.91517881772662e-9) * y + 8.35835975882941e-8) * y - 6.19785782240693e-7) *
				y + 3.95841149373135e-6) * y - 2.11366761402403e-5) * y + 9.00474771229507e-5) *
				y - 2.78777909813289e-4) * y + 5.26543779837487e-4;
		} else {	//x between 10 and 15
			y = x - 12.5;
			ldata->root[0] = (((((((((((y * 4.94869622744119e-17 + 8.0356880573916e-16) * y - 
				5.599125915431e-15) * y - 1.378685560217e-13) * y + 7.006511663249e-13) * y +
				1.30391406991118e-11) * y + 8.06987313467541e-11) * y - 5.20644072732933e-9) *
				y + 7.72794187755457e-8) * y - 1.61512612564194e-6) * y + 4.15083811185831e-5) *
				y - 7.87855975560199e-4) * y + .0114189319050009;
			ldata->root[1] = (((((((((((y * 4.89224285522336e-16 + 1.06390248099712e-14) * y - 
				5.446260182933e-14) * y - 1.613630106295e-12) * y + 3.910179118937e-12) * y +
				1.90712434258806e-10) * y + 8.78470199094761e-10) * y - 5.97332993206797e-8) *
				y + 9.25750831481589e-7) * y - 2.02362185197088e-5) * y + 4.92341968336776e-4) *
				y - .00868438439874703) * y + .115825965127958;
			ldata->root[2] = ((((((((((y * 6.12419396208408e-14 + 1.12328861406073e-13) * y - 
				9.051094103059e-12) * y - 4.781797525341e-11) * y + 1.660828868694e-9) * y +
				4.499058798868e-10) * y - 2.519549641933e-7) * y + 4.97744404018e-6) * y - 
				1.25858350034589e-4) * y + .00270279176970044) * y - .0399327850801083) * y +
				.433467200855434;
			ldata->root[3] = (((((((((((y * 4.63414725924048e-14 - 4.72757262693062e-14) * y - 
				1.001926833832e-11) * y + 6.074107718414e-11) * y + 1.576976911942e-9) * y -
				2.01186401974027e-8) * y - 1.84530195217118e-7) * y + 5.02333087806827e-6) * y +
				9.66961790843006e-6) * y - .00158522208889528) * y + .0280539673938339) * y -
				.278953904330072) * y + 1.82835655238235;
			ldata->weight[3] = (((((((((((((y * 2.90401781000996e-18 - 4.63389683098251e-17) * y 
				+ 6.274018198326e-16) * y - 8.936002188168e-15) * y + 1.194719074934e-13) * y -
				1.45501321259466e-12) * y + 1.64090830181013e-11) * y - 1.71987745310181e-10) *
				y + 1.63738403295718e-9) * y - 1.39237504892842e-8) * y + 1.06527318142151e-7) *
				y - 7.27634957230524e-7) * y + 4.12159381310339e-6) * y - 1.74648169719173e-5) *
				y + 8.50290130067818e-5;
			ldata->weight[2] = ((((((((((((y * -4.1956914545948e-17 + 5.94344180261644e-16) * y - 
				1.148797566469e-14) * y + 1.881303962576e-13) * y - 2.413554618391e-12) * y +
				3.372127423047e-11) * y - 4.933988617784e-10) * y + 6.116545396281e-9) * y - 
				6.69965691739299e-8) * y + 7.52380085447161e-7) * y - 8.08708393262321e-6) * y +
				6.88603417296672e-5) * y - 4.67067112993427e-4) * y + .00542313365864597;
			ldata->weight[1] = ((((((((((y * -6.22272689880615e-15 + 1.04126809657554e-13) * y - 
				6.842418230913e-13) * y + 1.576841731919e-11) * y - 4.203948834175e-10) * y +
				6.287255934781e-9) * y - 8.307159819228e-8) * y + 1.356478091922e-6) * y - 
				2.08065576105639e-5) * y + 2.5239673033234e-4) * y - .00294484050194539) * y +
				.0601396183129168;
			ldata->weight[0] = (((-.18784686463512 / x + .22991849164985) / x - .49893752514047)
				/ x - 2.1916512131607e-5) * exp(-x) + sqrt(pie4 / x) - ldata->weight[3] -
				ldata->weight[2] - ldata->weight[1];
		}
	} else {	// x > 15
		ldata->weight[0] = sqrt(pie4 / x);
		if (x <= 20.0) {	//x between 15 and 20
			y = x - 17.5;
			ldata->root[0] = (((((((((((y * 4.36701759531398e-17 - 1.12860600219889e-16) * y - 
				6.149849164164e-15) * y + 5.820231579541e-14) * y + 4.396602872143e-13) * y -
				1.24330365320172e-11) * y + 6.71083474044549e-11) * y + 2.43865205376067e-10) *
				y + 1.67559587099969e-8) * y - 9.32738632357572e-7) * y + 2.39030487004977e-5) *
				y - 4.68648206591515e-4) * y + .00834977776583956;
			ldata->root[1] = (((((((((((y * 4.98913142288158e-16 - 2.60732537093612e-16) * y - 
				7.775156445127e-14) * y + 5.766105220086e-13) * y + 6.4326967296e-12) * y -
				1.39571683725792e-10) * y + 5.95451479522191e-10) * y + 2.42471442836205e-9) *
				y + 2.4748571014312e-7) * y - 1.14710398652091e-5) * y + 2.71252453754519e-4) *
				y - .00496812745851408) * y + .082602060202678;
			ldata->root[2] = (((((((((((y * 1.91498302509009e-15 + 1.48840394311115e-14) * y - 
				4.316925145767e-13) * y + 1.186495793471e-12) * y + 4.615806713055e-11) * y -
				5.54336148667141e-10) * y + 3.48789978951367e-10) * y - 2.79188977451042e-9) *
				y + 2.09563208958551e-6) * y - 6.76512715080324e-5) * y + .00132129867629062) *
				y - .0205062147771513) * y + .288068671894324;
			ldata->root[3] = (((((((((((y * -5.43697691672942e-15 - 1.12483395714468e-13) * y + 
				2.826607936174e-12) * y - 1.26673449328e-11) * y - 4.258722866437e-10) * y +
				9.45486578503261e-9) * y - 5.86635622821309e-8) * y - 1.28835028104639e-6) * y + 
				4.41413815691885e-5) * y - 7.61738385590776e-4) * y + .0096609090298555) * y -
				.101410568057649) * y + .954714798156712;
			ldata->weight[3] = ((((((((((((y * -7.56882223582704e-19 + 7.53541779268175e-18) * y 
				- 1.157318032236e-16) * y + 2.411195002314e-15) * y - 3.601794386996e-14) * y +
				4.082150659615e-13) * y - 4.289542980767e-12) * y + 5.086829642731e-11) * y - 
				6.35435561050807e-10) * y + 6.82309323251123e-9) * y - 5.63374555753167e-8) *
				y + 3.57005361100431e-7) * y - 2.40050045173721e-6) * y + 4.94171300536397e-5;
			ldata->weight[2] = (((((((((((y * -5.54451040921657e-17 + 2.68748367250999e-16) * y + 
				1.349020069254e-14) * y - 2.507452792892e-13) * y + 1.944339743818e-12) * y -
				1.29816917658823e-11) * y + 3.49977768819641e-10) * y - 8.67270669346398e-9) *
				y + 1.31381116840118e-7) * y - 1.36790720600822e-6) * y + 1.1921069767316e-5) *
				y - 1.42181943986587e-4) * y + .00412615396191829;
			ldata->weight[1] = (((((((((((y * -1.865060577297e-16 + 1.16661114435809e-15) * y + 
				2.563712856363e-14) * y - 4.498350984631e-13) * y + 1.765194089338e-12) * y +
				9.04483676345625e-12) * y + 4.98930345609785e-10) * y - 2.11964170928181e-8) *
				y + 3.98295476005614e-7) * y - 5.49390160829409e-6) * y + 7.74065155353262e-5) *
				y - .00148201933009105) * y + .0497836392625268;
			ldata->weight[0] = ((.1962326414943 / x - .4969524146449) / x - 6.0156581186481e-5)
				* exp(-x) + ldata->weight[0] - ldata->weight[1] - ldata->weight[2] - 
				ldata->weight[3];
		} else if (x <= 35) {	//x between 20 and 35
			e = exp(-x);
			ldata->root[0] = ((((((x * -4.45711399441838e-5 + .00127267770241379) * x -
				.236954961381262) * x + 15.4330657903756) * x - 522.799159267808) * x +
				10595.1216669313) * x + (-2511772.35556236 / x + 872975.373557709) / 
				x - 129194.382386499) * e + r14 / (x - r14);
			ldata->root[1] = (((((x * -.0785617372254488 + 6.35653573484868) * x -
				338.29693876399) * x + 12512.0495802096) * x - 316847.570511637) * x +
				((-1024274661.27427 / x + 370104713.293016) / x - 58711900.5093822) / x +
				5386142.11391604) * e + r24 / (x - r24);
			ldata->root[2] = (((((x * -.237900485051067 + 18.4122184400896) * x -
				1002.00731304146) * x + 37515.1841595736) * x - 950626.66339013) * x +
				((-2881390146.51985 / x + 1066259150.44526) / x - 172465289.687396) / x + 
				16041939.0230055) * e + r34 / (x - r34);
			ldata->root[3] = ((((((x * -6.00691586407385e-4 - .364479545338439) * x +
				15.7496131755179) * x - 654.944248734901) * x + 17083.0039597097) * x -
				290517.939780207) * x + (34905969.8304732 / x - 16494452.2586065) / x +
				2968179.40164703) * e + r44 / (x - r44);
			if (x <= 25.) {
				ldata->weight[3] = (((((((x * 2.33766206773151e-7 - 3.81542906607063e-5) * 
					x + .00351416601267) * x - .166538571864728) * x + 4.80006136831847) *
					x - 87.3165934223603) * x + 977.683627474638) * x + 16600.094511764 / 
					x - 6144.79071209961) * e + w44 * ldata->weight[0];
			} else {
				ldata->weight[3] = ((((((x * 5.74245945342286e-6 - 7.58735928102351e-5) * 
					x + 2.35072857922892e-4) * x - .00378812134013125) * x + .309871652785805) *
					x - 7.11108633061306) * x + 55.5297573149528) * e + w44 * ldata->weight[0];
			}
			ldata->weight[2] = ((((((x * 2.36392855180768e-4 - .00916785337967013) * x +
				.462186525041313) * x - 19.694378600654) * x + 499.169195295559) * x -
				6214.1984584509) * x + ((52144505.3212414 / x - 13411346.4389309) / x +
				1136732.98305631) / x - 2815.01182042707) * e + w34 * ldata->weight[0];
			ldata->weight[1] = ((((((x * 7.29841848989391e-4 - .0353899555749875) * x +
				2.07797425718513) * x - 100.464709786287) * x + 3152.06108877819) * x -
				62705.4715090012) * x + (15472124.6264919 / x - 5260743.91316381) / x +
				767135.400969617) * e + w24 * ldata->weight[0];
			ldata->weight[0] = ((.1962326414943 / x - .4969524146449) / x - 6.0156581186481e-5)
				* e + ldata->weight[0] - ldata->weight[1] - ldata->weight[2] - ldata->weight[3];
		} else if (x <= 53.0) {	//x between 35 and 53
			double	x2;
			x2 = x * x;
			e = exp(-x) * (x2 * x2);
			ldata->root[3] = ((x * -.00219135070169653 - .119108256987623) * x -
				.750238795695573) * e + r44 / (x - r44);
			ldata->root[2] = ((x * -9.65842534508637e-4 - .0449822013469279) * x 
				+ .608784033347757) * e + r34 / (x - r34);
			ldata->root[1] = ((x * -3.62569791162153e-4 - .00909231717268466) * x 
				+ .184336760556262) * e + r24 / (x - r24);
			ldata->root[0] = ((x * -4.075575259146e-5 - 6.88846864931685e-4) * x 
				+ .0174725309199384) * e + r14 / (x - r14);
			ldata->weight[3] = ((x * 5.7663198200099e-6 - 7.8918728380489e-5) * x + 
				3.28297971853126e-4) * e + w44 * ldata->weight[0];
			ldata->weight[2] = ((x * 2.0829496985723e-4 - .00377489954837361) * x + 
				.0209857151617436) * e + w34 * ldata->weight[0];
			ldata->weight[1] = ((x * 6.16374517326469e-4 - .0126711744680092) * x + 
				.0814504890732155) * e + w24 * ldata->weight[0];
			ldata->weight[0] -= ldata->weight[1] + ldata->weight[2] + ldata->weight[3];
		} else {	// x > 53
			ldata->root[0] = r14 / (x - r14);
			ldata->root[1] = r24 / (x - r24);
			ldata->root[2] = r34 / (x - r34);
			ldata->root[3] = r44 / (x - r44);
			ldata->weight[3] = w44 * ldata->weight[0];
			ldata->weight[2] = w34 * ldata->weight[0];
			ldata->weight[1] = w24 * ldata->weight[0];
			ldata->weight[0] -= ldata->weight[1] + ldata->weight[2] + ldata->weight[3];
		}
	}
}

#define r15		.117581320211778
#define r25		1.0745620124369
#define w25		.270967405960535
#define r35		3.08593744371754
#define w35		.0382231610015404
#define r45		6.41472973366203
#define w45		.00151614186862443
#define r55		11.8071894899717
#define w55		8.62130526143657e-6

void Root5(Root * ldata) {
	double e, y;

/*          *****   VERSION  FEBRUARY 27,1975   ***** */

	double x = ldata->x;
	if (x <= 15.) {
		if (x <= 3.e-7) {	//x is almost zero
			ldata->root[0] = .0226659266316985 - x * .00215865967920897;
			ldata->root[1] = .231271692140903 - x * .0220258754389745;
			ldata->root[2] = .857346024118836 - x * .0816520023025515;
			ldata->root[3] = 2.97353038120346 - x * .283193369647137;
			ldata->root[4] = 18.4151859759051 - x * 1.75382723579439;
			ldata->weight[0] = .295524224714752 - x * .0196867576909777;
			ldata->weight[1] = .269266719309995 - x * .0561737590184721;
			ldata->weight[2] = .219086362515981 - x * .0971152726793658;
			ldata->weight[3] = .14945134915058 - x * .102979262193565;
			ldata->weight[4] = .0666713443086877 - x * .0573782817488315;
		} else if (x <= 1.0) {	//x between zero and 1
			ldata->root[0] = ((((((x * -4.46679165328413e-11 + 1.21879111988031e-9) * x -
				2.62975022612104e-8) * x + 5.15106194905897e-7) * x - 9.27933625824749e-6) *
				x + 1.51794097682482e-4) * x - .00215865967920301) * x + .0226659266316985;
			ldata->root[1] = ((((((x * 1.93117331714174e-10 - 4.57267589660699e-9) * x +
				2.48339908218932e-8) * x + 1.50716729438474e-6) * x - 6.07268757707381e-5) *
				x + .00137506939145643) * x - .0220258754419939) * x + .231271692140905;
			ldata->root[2] = (((((x * 4.84989776180094e-9 + 1.31538893944284e-7) * x -
				2.766753852879e-6) * x - 7.651163510626e-5) * x + .004033058545972) * x -
				.0816520022916145) * x + .857346024118779;
			ldata->root[3] = ((((x * -2.48581772214623e-7 - 4.34482635782585e-6) * x -
				7.4601825798763e-7) * x + .0101210776517279) * x - .283193369640005) * x +
				2.97353038120345;
			ldata->root[4] = (((((x * -8.92432153868554e-9 + 1.77288899268988e-8) * x +
				3.040754680666e-6) * x + 1.058229325071e-4) * x + .04596379534985) * x -
				1.75382723579114) * x + 18.4151859759049;
			ldata->weight[0] = ((((((x * -2.03822632771791e-9 + 3.8911022913381e-8) * x -
				5.84914787904823e-7) * x + 8.30316168666696e-6) * x - 1.13218402310546e-4) *
				x + .0014912888858679) * x - .0196867576904816) * x + .295524224714749;
			ldata->weight[1] = (((((((x * 8.6284811839757e-9 - 1.38975551148989e-7) * x +
				1.602894068228e-6) * x - 1.646364300836e-5) * x + 1.538445806778e-4) * x -
				.00128848868034502) * x + .00938866933338584) * x - .0561737590178812) * 
				x + .269266719309991;
			ldata->weight[2] = ((((((((x * -9.41953204205665e-9 + 1.47452251067755e-7) * x -
				1.57456991199322e-6) * x + 1.45098401798393e-5) * x - 1.18858834181513e-4) *
				x + 8.5369767598421e-4) * x - .00522877807397165) * x + .0260854524809786) *
				x - .0971152726809059) * x + .219086362515979;
			ldata->weight[3] = ((((((((x * -3.84961617022042e-8 + 5.6659539654447e-7) * x -
				5.52351805403748e-6) * x + 4.53160377546073e-5) * x - 3.22542784865557e-4) *
				x + .00195682017370967) * x - .00977232537679229) * x + .0379455945268632) *
				x - .102979262192227) * x + .149451349150573;
			ldata->weight[4] = (((((((((x * 4.0959481252143e-9 - 6.47097874264417e-8) * x +
				6.743541482689e-7) * x - 5.917993920224e-6) * x + 4.531969237381e-5) * x -
				2.99102856679638e-4) * x + .00165695765202643) * x - .00740671222520653) * 
				x + .0250889946832192) * x - .0573782817487958) * x + .0666713443086877;
		} else if (x <= 5.) {	//x between 1 and 5
			y = x - 3.;
			ldata->root[0] = ((((((((y * -2.58163897135138e-14 + 8.14127461488273e-13) * y - 
				2.11414838976129e-11) * y + 5.09822003260014e-10) * y - 1.16002134438663e-8) *
				y + 2.4681069441454e-7) * y - 4.92556826124502e-6) * y + 9.02580687971053e-5) *
				y - .00145190025120726) * y + .0173416786387475;
			ldata->root[1] = (((((((((y * 1.04525287289788e-14 + 5.44611782010773e-14) * y - 
				4.831059411392e-12) * y + 1.136643908832e-10) * y - 1.104373076913e-9) * y -
				2.35346740649916e-8) * y + 1.43772622028764e-6) * y - 4.23405023015273e-5) *
				y + 9.12034574793379e-4) * y - .0152479441718739) * y + .176055265928744;
			ldata->root[2] = (((((((((y * -6.89693150857911e-14 + 5.92064260918861e-13) * y + 
				1.847170956043e-11) * y - 3.390752744265e-10) * y - 2.995532064116e-9) * y +
				1.57456141058535e-7) * y - 3.95859409711346e-7) * y - 9.58924580919747e-5) *
				y + .00323551502557785) * y - .0597587007636479) * y + .646432853383057;
			ldata->root[3] = ((((((((y * -3.61293809667763e-12 - 2.70803518291085e-11) * y + 
				8.83758848468769e-10) * y + 1.59166632851267e-8) * y - 1.32581997983422e-7) *
				y - 7.60223407443995e-6) * y - 7.41019244900952e-5) * y + .00981432631743423) *
				y - .223055570487771) * y + 2.21460798080643;
			ldata->root[4] = (((((((((y * 7.12332088345321e-13 + 3.16578501501894e-12) * y - 
				8.776668218053e-11) * y - 2.342817613343e-9) * y - 3.496962018025e-8) * y -
				3.03172870136802e-7) * y + 1.50511293969805e-6) * y + 1.37704919387696e-4) *
				y + .0470723869619745) * y - 1.47486623003693) * y + 13.5704792175847;
			ldata->weight[0] = (((((((((y * 1.04348658616398e-13 - 1.94147461891055e-12) * y + 
				3.485512360993e-11) * y - 6.277497362235e-10) * y + 1.100758247388e-8) * y -
				1.88329804969573e-7) * y + 3.12338120839468e-6) * y - 5.04404167403568e-5) *
				y + 8.00338056610995e-4) * y - .0130892406559521) * y + .247383140241103;
			ldata->weight[1] = (((((((((((y * 3.23496149760478e-14 - 5.24314473469311e-13) * y + 
				7.743219385056e-12) * y - 1.146022750992e-10) * y + 1.615238462197e-9) * y -
				2.15479017572233e-8) * y + 2.70933462557631e-7) * y - 3.18750295288531e-6) * y + 
				3.47425221210099e-5) * y - 3.45558237388223e-4) * y + .00305779768191621) * y -
				.0229118251223003) * y + .159834227924213;
			ldata->weight[2] = ((((((((((((y * -3.42790561802876e-14 + 5.26475736681542e-13) * y 
				- 7.184330797139e-12) * y + 9.763932908544e-11) * y - 1.244014559219e-9) * y +
				1.472744068942e-8) * y - 1.611749975234e-7) * y + 1.616487851917e-6) * y - 
				1.46852359124154e-5) * y + 1.18900349101069e-4) * y - 8.37562373221756e-4) *
				y + .00493752683045845) * y - .0225514728915673) * y + .0695211812453929;
			ldata->weight[3] = (((((((((((((y * 1.04072340345039e-14 - 1.60808044529211e-13) * y 
				+ 2.183534866798e-12) * y - 2.939403008391e-11) * y + 3.679254029085e-10) * y -
				4.23775673047899e-9) * y + 4.46559231067006e-8) * y - 4.26488836563267e-7) * y + 
				3.64721335274973e-6) * y - 2.74868382777722e-5) * y + 1.78586118867488e-4) *
				y - 9.68428981886534e-4) * y + .00416002324339929) * y - .0128290192663141) *
				y + .0222353727685016;
			ldata->weight[4] = ((((((((((((((y * -8.16770412525963e-16 + 1.31376515047977e-14) * 
				y - 1.856950818865e-13) * y + 2.596836515749e-12) * y - 3.372639523006e-11) *
				y + 4.025371849467e-10) * y - 4.389453269417e-9) * y + 4.332753856271e-8) * y - 
				3.82673275931962e-7) * y + 2.98006900751543e-6) * y - 2.00718990300052e-5) *
				y + 1.13876001386361e-4) * y - 5.23627942443563e-4) * y + .00183524565118203) *
				y - .00437785737450783) * y + .00536963805223095;
		} else if (x <= 10.) {	//x between 5 and 10
			y = x - 7.5;
			ldata->root[0] = ((((((((y * -1.13825201010775e-14 + 1.89737681670375e-13) * y - 
				4.81561201185876e-12) * y + 1.56666512163407e-10) * y - 3.73782213255083e-9) *
				y + 9.15858355075147e-8) * y - 2.13775073585629e-6) * y + 4.56547356365536e-5) *
				y - 8.6800390932374e-4) * y + .0122703754069176;
			ldata->root[1] = (((((((((y * -3.67160504428358e-15 + 1.27876280158297e-14) * y - 
				1.296476623788e-12) * y + 1.477175434354e-11) * y + 5.464102147892e-10) * y -
				2.42538340602723e-8) * y + 8.20460740637617e-7) * y - 2.20379304598661e-5) * y + 
				4.90295372978785e-4) * y - .00914294111576119) * y + .12259040340369;
			ldata->root[2] = (((((((((y * 1.39017367502123e-14 - 6.9639138542689e-13) * y + 
				1.176946020731e-12) * y + 1.725627235645e-10) * y - 3.6863838563e-9) * y +
				2.87495324207095e-8) * y + 1.71307311000282e-6) * y - 7.94273603184629e-5) *
				y + .00200938064965897) * y - .0363329491677178) * y + .434393683888443;
			ldata->root[3] = ((((((((((y * -1.27815158195209e-14 + 1.99910415869821e-14) * y + 
				3.753542914426e-12) * y - 2.708018219579e-11) * y - 1.190574776587e-9) * y +
				1.106696436509e-8) * y + 3.954955671326e-7) * y - 4.398596059588e-6) * y - 
				2.01087998907735e-4) * y + .00789092425542937) * y - .142056749162695) * y +
				1.39964149420683;
			ldata->root[4] = ((((((((((y * -1.19442341030461e-13 - 2.34074833275956e-12) * y + 
				6.861649627426e-12) * y + 6.082671496226e-10) * y + 5.38116010542e-9) * y -
				6.2532971387e-8) * y - 2.13596683505e-6) * y - 2.373394341886e-5) * y +
				2.88711171412814e-6) * y + .0485221195290753) * y - 1.04346091985269) * y +
				7.89901551676692;
			ldata->weight[0] = (((((((((y * 7.95526040108997e-15 - 2.48593096128045e-13) * y + 
				4.76124620872e-12) * y - 9.535763686605e-11) * y + 2.225273630974e-9) * y -
				4.49796778054865e-8) * y + 9.17812870287386e-7) * y - 1.86764236490502e-5) *
				y + 3.76807779068053e-4) * y - .00810456360143408) * y + .201097936411496;
			ldata->weight[1] = (((((((((((y * 1.25678686624734e-15 - 2.34266248891173e-14) * y + 
				3.973252415832e-13) * y - 6.830539401049e-12) * y + 1.140771033372e-10) * y -
				1.82546185762009e-9) * y + 2.77209637550134e-8) * y - 4.01726946190383e-7) * y + 
				5.48227244014763e-6) * y - 6.95676245982121e-5) * y + 8.05193921815776e-4) * y -
				.00815528438784469) * y + .0971769901268114;
			ldata->weight[2] = ((((((((((((y * -8.20929494859896e-16 + 1.37356038393016e-14) * y 
				- 2.02286306522e-13) * y + 3.058055403795e-12) * y - 4.387890955243e-11) * y +
				5.923946274445e-10) * y - 7.503659964159e-9) * y + 8.851599803902e-8) * y - 
				9.65561998415038e-7) * y + 9.60884622778092e-6) * y - 8.56551787594404e-5) * y +
				6.66057194311179e-4) * y - .00417753183902198) * y + .0225443826852447;
			ldata->weight[3] = ((((((((((((((y * -1.0876461248879e-17 + 1.85299909689937e-16) * y 
				- 2.730195628655e-15) * y + 4.127368817265e-14) * y - 5.881379088074e-13) * y +
				7.805245193391e-12) * y - 9.632707991704e-11) * y + 1.099047050624e-9) * y - 
				1.15042731790748e-8) * y + 1.09415155268932e-7) * y - 9.33687124875935e-7) * y +
				7.02338477986218e-6) * y - 4.53759748787756e-5) * y + 2.41722511389146e-4) * y - 
				9.75935943447037e-4) * y + .00257520532789644;
			ldata->weight[4] = (((((((((((((((y * 7.28996979748849e-19 - 1.26518146195173e-17) * 
				y + 1.886145834486e-16) * y - 2.876728287383e-15) * y + 4.114588668138e-14) *
				y - 5.44436631413933e-13) * y + 6.64976446790959e-12) * y - 7.4456006997494e-11) *
				y + 7.57553198166848e-10) * y - 6.92956101109829e-9) * y + 5.62222859033624e-8) *
				y - 3.97500114084351e-7) * y + 2.3903912613814e-6) * y - 1.18023950002105e-5) *
				y + 4.52254031046244e-5) * y - 1.2111378215037e-4) * y + 1.75013126731224e-4;
		} else {	//x between 10 and 15
			y = x - 12.5;
			ldata->root[0] = ((((((((((y * -4.16387977337393e-17 + 7.2087299737386e-16) * y + 
				1.395993802064e-14) * y + 3.660484641252e-14) * y - 4.154857548139e-12) * y +
				2.301379846544e-11) * y - 1.033307012866e-9) * y + 3.997777641049e-8) * y - 
				9.35118186333939e-7) * y + 2.38589932752937e-5) * y - 5.35185183652937e-4) *
				y + .00885218988709735;
			ldata->root[1] = ((((((((((y * -4.56279214732217e-16 + 6.24941647247927e-15) * y + 
				1.737896339191e-13) * y + 8.964205979517e-14) * y - 3.538906780633e-11) * y +
				9.561341254948e-11) * y - 9.77283189131e-9) * y + 4.24034019462e-7) * y - 
				1.02384302866534e-5) * y + 2.57987709704822e-4) * y - .00554735977651677) *
				y + .0868245143991948;
			ldata->root[2] = ((((((((((y * -2.52879337929239e-15 + 2.13925810087833e-14) * y + 
				7.884307667104e-13) * y - 9.02339815951e-13) * y - 5.814101544957e-11) * y -
				1.333480437968e-9) * y - 2.217064940373e-8) * y + 1.643290788086e-6) * y - 
				4.39602147345028e-5) * y + .00108648982748911) * y - .0213014521653498) *
				y + .294150684465425;
			ldata->root[3] = ((((((((((y * -6.42391438038888e-15 + 5.37848223438815e-15) * y + 
				8.960828117859e-13) * y + 5.214153461337e-11) * y - 1.106601744067e-10) * y -
				2.007890743962e-8) * y + 1.543764346501e-7) * y + 4.520749076914e-6) * y - 
				1.88893338587047e-4) * y + .00473264487389288) * y - .0791197893350253) * y +
				.860057928514554;
			ldata->root[4] = (((((((((((y * -2.24366166957225e-14 + 4.87224967526081e-14) * y + 
				5.587369053655e-12) * y - 3.045253104617e-12) * y - 1.22398388308e-9) * y -
				2.05603889396319e-9) * y + 2.58604071603561e-7) * y + 1.34240904266268e-6) * y - 
				5.72877569731162e-5) * y - 9.56275105032191e-4) * y + .0423367010370921) * y -
				.576800927133412) * y + 3.87328263873381;
			ldata->weight[0] = (((((((((y * 8.98007931950169e-15 + 7.25673623859497e-14) * y + 
				5.851494250405e-14) * y - 4.234204823846e-11) * y + 3.911507312679e-10) * y -
				9.65094802088511e-9) * y + 3.42197444235714e-7) * y - 7.51821178144509e-6) * y + 
				1.94218051498662e-4) * y - .00538533819142287) * y + .168122596736809;
			ldata->weight[1] = ((((((((((y * -1.05490525395105e-15 + 1.96855386549388e-14) * y - 
				5.500330153548e-13) * y + 1.003849567976e-11) * y - 1.720997242621e-10) * y +
				3.533277061402e-9) * y - 6.389171736029e-8) * y + 1.046236652393e-6) * y - 
				1.73148206795827e-5) * y + 2.57820531617185e-4) * y - .0034618826533835) *
				y + .0703302497508176;
			ldata->weight[2] = (((((((((((y * 3.60020423754545e-16 - 6.24245825017148e-15) * y + 
				9.945311467434e-14) * y - 1.749051512721e-12) * y + 2.768503957853e-11) * y -
				4.08688551136506e-10) * y + 6.0418906330361e-9) * y - 8.23540111024147e-8) * y + 
				1.01503783870262e-6) * y - 1.20490761741576e-5) * y + 1.26928442448148e-4) *
				y - .00105539461930597) * y + .0115543698537013;
			ldata->weight[3] = (((((((((((((y * 2.51163533058925e-18 - 4.31723745510697e-17) * y 
				+ 6.557620865832e-16) * y - 1.016528519495e-14) * y + 1.491302084832e-13) * y -
				2.06638666222265e-12) * y + 2.67958697789258e-11) * y - 3.23322654638336e-10) *
				y + 3.63722952167779e-9) * y - 3.75484943783021e-8) * y + 3.49164261987184e-7) *
				y - 2.92658670674908e-6) * y + 2.12937256719543e-5) * y - 1.19434130620929e-4) *
				y + 6.45524336158384e-4;
			ldata->weight[4] = ((((((((((((((y * -1.29043630202811e-19 + 2.16234952241296e-18) * 
				y - 3.107631557965e-17) * y + 4.570804313173e-16) * y - 6.301348858104e-15) *
				y + 8.031304476153e-14) * y - 9.446196472547e-13) * y + 1.018245804339e-11) *
				y - 9.96995451348129e-11) * y + 8.77489010276305e-10) * y - 6.84655877575364e-9)
				* y + 4.64460857084983e-8) * y - 2.66924538268397e-7) * y + 1.24621276265907e-6)
				* y - 4.30868944351523e-6) * y + 9.94307982432868e-6;
		}
	} else {
		if (x <= 20.) {	//x between 15 and 20
			y = x - 17.5;
			ldata->root[0] = ((((((((((y * 1.9187576454574e-16 + 7.8357401095707e-16) * y - 
				3.260875931644e-14) * y - 1.186752035569e-13) * y + 4.275180095653e-12) * y +
				3.357056136731e-11) * y - 1.123776903884e-9) * y + 1.231203269887e-8) * y - 
				3.99851421361031e-7) * y + 1.45418822817771e-5) * y - 3.49912254976317e-4) *
				y + .00667768703938812;
			ldata->root[1] = ((((((((((y * 2.02778478673555e-15 + 1.01640716785099e-14) * y - 
				3.385363492036e-13) * y - 1.615655871159e-12) * y + 4.527419140333e-11) * y +
				3.853670706486e-10) * y - 1.184607130107e-8) * y + 1.347873288827e-7) * y - 
				4.47788241748377e-6) * y + 1.54942754358273e-4) * y - .00355524254280266) *
				y + .0644912219301603;
			ldata->root[2] = ((((((((((y * 7.79850771456444e-15 + 6.00464406395001e-14) * y - 
				1.249779730869e-12) * y - 1.020720636353e-11) * y + 1.814709816693e-10) * y +
				1.766397336977e-9) * y - 4.60355944901e-8) * y + 5.863956443581e-7) * y - 
				2.03797212506691e-5) * y + 6.31405161185185e-4) * y - .0130102750145071) *
				y + .210244289044705;
			ldata->root[3] = (((((((((((y * -2.92397030777912e-15 + 1.94152129078465e-14) * y + 
				4.85944766585e-13) * y - 3.217227223463e-12) * y - 7.484522135512e-11) * y +
				7.19101516047753e-10) * y + 6.88409355245582e-9) * y - 1.44374545515769e-7) *
				y + 2.74941013315834e-6) * y - 1.02790452049013e-4) * y + .00259924221372643) *
				y - .0435712368303551) * y + .562170709585029;
			ldata->root[4] = (((((((((((y * 1.1797612684006e-14 + 1.24156229350669e-13) * y - 
				3.89274162228e-12) * y - 7.755793199043e-12) * y + 9.492190032313e-10) * y -
				4.98680128123353e-9) * y - 1.81502268782664e-7) * y + 2.69463269394888e-6) *
				y + 2.5003215442164e-5) * y - .00133684303917681) * y + .0229121951862538) *
				y - .245653725061323) * y + 1.89999883453047;
			ldata->weight[0] = ((((((((((y * 1.74841995087592e-15 - 6.95671892641256e-16) * y - 
				3.000659497257e-13) * y + 2.021279817961e-13) * y + 3.8535969354e-11) * y +
				1.461418533652e-10) * y - 1.014517563435e-8) * y + 1.132736008979e-7) * y - 
				2.86605475073259e-6) * y + 1.21958354908768e-4) * y - .00386293751153466) *
				y + .145298342081522;
			ldata->weight[1] = ((((((((((y * -1.11199320525573e-15 + 1.85007587796671e-15) * y + 
				1.220613939709e-13) * y + 1.275068098526e-12) * y - 5.341838883262e-11) * y +
				6.161037256669e-10) * y - 1.00914787975e-8) * y + 2.907862965346e-7) * y - 
				6.12300038720919e-6) * y + 1.00104454489518e-4) * y - .00180677298502757) *
				y + .057800991453663;
			ldata->weight[2] = ((((((((((y * -9.49816486853687e-16 + 6.67922080354234e-15) * y + 
				2.606163540537e-15) * y + 1.98379995015e-12) * y - 5.400548574357e-11) * y +
				6.638043374114e-10) * y - 8.799518866802e-9) * y + 1.791418482685e-7) * y - 
				2.96075397351101e-6) * y + 3.38028206156144e-5) * y - 3.58426847857878e-4) *
				y + .00839213709428516;
			ldata->weight[3] = (((((((((((y * 1.3382997106018e-17 - 3.4484187784414e-16) * y + 
				4.745009557656e-15) * y - 6.033814209875e-14) * y + 1.049256040808e-12) * y -
				1.70859789556117e-11) * y + 2.15219425727959e-10) * y - 2.52746574206884e-9) *
				y + 3.2776171442296e-8) * y - 3.90387662925193e-7) * y + 3.4634020459387e-6) *
				y - 2.43236345136782e-5) * y + 3.54846978585226e-4;
			ldata->weight[4] = (((((((((((((y * 2.69412277020887e-20 - 4.24837886165685e-19) * y 
				+ 6.030500065438e-18) * y - 9.069722758289e-17) * y + 1.246599177672e-15) * y -
				1.56872999797549e-14) * y + 1.87305099552692e-13) * y - 2.09498886675861e-12) *
				y + 2.11630022068394e-11) * y - 1.92566242323525e-10) * y + 1.62012436344069e-9)
				* y - 1.23621614171556e-8) * y + 7.72165684563049e-8) * y - 3.59858901591047e-7)
				* y + 2.43682618601e-6;
		} else if (x <= 25.) {	//x between 20 and 25
			y = x - 22.5;
			ldata->root[0] = (((((((((y * -1.13927848238726e-15 + 7.39404133595713e-15) * y + 
				1.445982921243e-13) * y - 2.676703245252e-12) * y + 5.823521627177e-12) * y +
				2.17264723874381e-10) * y + 3.56242145897468e-9) * y - 3.03763737404491e-7) *
				y + 9.46859114120901e-6) * y - 2.30896753853196e-4) * y + .00524663913001114;
			ldata->root[1] = ((((((((((y * 2.89872355524581e-16 - 1.22296292045864e-14) * y + 
				6.1840650972e-14) * y + 1.64984659123e-12) * y - 2.729713905266e-11) * y +
				3.70991379065e-11) * y + 2.216486288382e-9) * y + 4.616160236414e-8) * y - 
				3.32380270861364e-6) * y + 9.84635072633776e-5) * y - .00230092118015697) *
				y + .0500845183695073;
			ldata->root[2] = ((((((((((y * 1.97068646590923e-15 - 4.894192706268e-14) * y + 
				1.136466605916e-13) * y + 7.546203883874e-12) * y - 9.635646767455e-11) * y -
				8.295965491209e-11) * y + 7.534109114453e-9) * y + 2.699970652707e-7) * y - 
				1.42982334217081e-5) * y + 3.78290946669264e-4) * y - .00803133015084373) *
				y + .158689469640791;
			ldata->root[3] = ((((((((((y * 1.33642069941389e-14 - 1.55850612605745e-13) * y - 
				7.522712577474e-13) * y + 3.209520801187e-11) * y - 2.075594313618e-10) * y -
				2.070575894402e-9) * y + 7.323046997451e-9) * y + 1.851491550417e-6) * y - 
				6.37524802411383e-5) * y + .00136795464918785) * y - .0242051126993146) *
				y + .397847167557815;
			ldata->root[4] = ((((((((((y * -6.07053986130526e-14 + 1.04447493138843e-12) * y - 
				4.286617818951e-13) * y - 2.632066100073e-10) * y + 4.804518986559e-9) * y -
				1.835675889421e-8) * y - 1.068175391334e-6) * y + 3.292234974141e-5) * y - 
				5.94805357558251e-4) * y + .00829382168612791) * y - .0993122509049447) *
				y + 1.09857804755042;
			ldata->weight[0] = (((((((((y * -9.10338640266542e-15 + 1.00438927627833e-13) * y + 
				7.817349237071e-13) * y - 2.547619474232e-11) * y + 1.479321506529e-10) * y +
				1.52314028857627e-9) * y + 9.20072040917242e-9) * y - 2.19427111221848e-6) * y + 
				8.65797782880311e-5) * y - .00282718629312875) * y + .128718310443295;
			ldata->weight[1] = (((((((((y * 5.5238092761876e-15 - 6.43424400204124e-14) * y - 
				2.358734508092e-13) * y + 8.261326648131e-12) * y + 9.229645304956e-11) * y -
				5.68108973828949e-9) * y + 1.22477891136278e-7) * y - 2.11919643127927e-6) * y +
				4.23605032368922e-5) * y - .00114423444576221) * y + .0506607252890186;
			ldata->weight[2] = (((((((((y * 3.99457454087556e-15 - 5.11826702824182e-14) * y - 
				4.157593182747e-14) * y + 4.214670817758e-12) * y + 6.705582751532e-11) * y -
				3.36086411698418e-9) * y + 6.07453633298986e-8) * y - 7.40736211041247e-7) *
				y + 8.84176371665149e-6) * y - 1.72559275066834e-4) * y + .00716639814253567;
			ldata->weight[3] = (((((((((((y * -2.14649508112234e-18 - 2.45525846412281e-18) * y + 
				6.126212599772e-16) * y - 8.526651626939e-15) * y + 4.826636065733e-14) * y -
				3.3955416364974e-13) * y + 1.67070784862985e-11) * y - 4.42671979311163e-10) *
				y + 6.773680559084e-9) * y - 7.03520999708859e-8) * y + 6.04993294708874e-7) *
				y - 7.80555094280483e-6) * y + 2.85954806605017e-4;
			ldata->weight[4] = ((((((((((((y * -5.63938733073804e-21 + 6.92182516324628e-20) * y 
				- 1.586937691507e-18) * y + 3.357639744582e-17) * y - 4.810285046442e-16) * y +
				5.386312669975e-15) * y - 6.117895297439e-14) * y + 8.441808227634e-13) * y - 
				1.18527596836592e-11) * y + 1.36296870441445e-10) * y - 1.17842611094141e-9) *
				y + 7.80430641995926e-9) * y - 5.9776741740054e-8) * y + 1.65186146094969e-6;
		} else {
			ldata->weight[0] = sqrt(pie4 / x);
			if (x <= 40.) {	//x between 25 and 40
				e = exp(-x);
				ldata->root[0] = ((((((((x * -1.73363958895356e-6 + 1.19921331441483e-4) * 
					x - .0159437614121125) * x + 1.13467897349442) * x - 44.7216460864586) *
					x + 1062.51216612604) * x - 15207.3917378512) * x + 120662.887111273) * 
					x - 407186.366852475) * e + r15 / (x - r15);
				ldata->root[1] = ((((((((x * -1.6010254262171e-5 + .00110331262112395) * 
					x - .150043662589017) * x + 10.5563640866077) * x - 410.468817024806) *
					x + 9626.04416506819) * x - 135888.06983827) * x + 1061075.7703834) * 
					x - 3511907.92816119) * e + r25 / (x - r25);
				ldata->root[2] = ((((((((x * -4.48880032128422e-5 + .00269025112122177) * 
					x - .401048115525954) * x + 27.8360021977405) * x - 1048.91729356965) *
					x + 23698.5942687423) * x - 319504.627257548) * x + 2348796.93563358) * 
					x - 7163415.68174085) * e + r35 / (x - r35);
				ldata->root[3] = ((((((((x * -6.38526371092582e-5 - .00229263585792626) * 
					x - .0765735935499627) * x + 9.12692349152792) * x - 232.077034386717) *
					x + 281.839578728845) * x + 95952.9683876419) * x - 1776389.56809518) * 
					x + 10248975.964541) * e + r45 / (x - r45);
				ldata->root[4] = ((((((((x * -3.59049364231569e-5 - .0225963977930044) * 
					x + 1.12594870794668) * x - 45.6752462103909) * x + 1058.04526830637) *
					x - 11600.3199605875) * x - 40729.7627297272) * x + 2222155.28319857) * 
					x - 16119645.5032613) * e + r55 / (x - r55);
				ldata->weight[4] = (((((((((x * -4.6110090613397e-10 + 1.43069932644286e-7) * 
					x - 1.6396091543108e-5) * x + .00115791154612838) * x - .0530573476742071) *
					x + 1.61156533367153) * x - 32.3248143316007) * x + 412.007318109157) * 
					x - 3022.60070158372) * x + 9715.75094154768) * e + w55 * ldata->weight[0];
				ldata->weight[3] = (((((((((x * -2.4079943580995e-8 + 8.12621667601546e-6) * 
					x - 9.04491430884113e-4) * x + .0637686375770059) * x - 2.96135703135647) *
					x + 91.514235699633) * x - 1869.71865249111) * x + 24294.5528916947) * 
					x - 181852.473229081) * x + 596854.758661427) * e + w45 * ldata->weight[0];
				ldata->weight[2] = ((((((((x * 1.83574464457207e-5 - .00154837969489927) * 
					x + .118520453711586) * x - 6.69649981309161) * x + 244.789386487321) *
					x - 5688.32664556359) * x + 81450.7604229357) * x - 655181.056671474) * 
					x + 2264108.96607237) * e + w35 * ldata->weight[0];
				ldata->weight[1] = ((((((((x * 2.7777834587065e-5 - .0022283501765589) * 
					x + .161077633475573) * x - 8.96743743396132) * x + 328.062687293374) *
					x - 7657.22701219557) * x + 110255.055017664) * x - 892528.122219324) * 
					x + 3106386.27744347) * e + w25 * ldata->weight[0];
				ldata->weight[0] = ldata->weight[0] - e * .01962 - ldata->weight[1] -
					ldata->weight[2] - ldata->weight[3] - ldata->weight[4];
			} else if (x <= 59.) {	//x between 40 and 59
				double	x3;
				x3 = x * x * x;
				e = x3 * exp(-x);
				ldata->root[0] = (((x * -.0243758528330205 + 2.07301567989771) * x - 
					64.5964225381113) * x + 714.16008865547) * e + r15 / (x - r15);
				ldata->root[1] = (((x * -.228861955413636 + 19.3190784733691) * x - 
					599.774730340912) * x + 6618.44165304871) * e + r25 / (x - r25);
				ldata->root[2] = (((x * -.695053039285586 + 57.6874090316016) * x - 
					1777.0414322552) * x + 19536.6082947811) * e + r35 / (x - r35);
				ldata->root[3] = (((x * -1.58072809087018 + 127.050801091948) * x - 
					3866.8735091428) * x + 42302.482812142) * e + r45 / (x - r45);
				ldata->root[4] = (((x * -3.33963830405396 + 251.830424600204) * x - 
					7577.28527654961) * x + 82196.681659569) * e + r55 / (x - r55);
				e = x3 * e;
				ldata->weight[4] = ((x * 1.35482430510942e-8 - 3.27722199212781e-7) * x 
					+ 2.41522703684296e-6) * e + w55 * ldata->weight[0];
				ldata->weight[3] = ((x * 1.23464092261605e-6 - 3.5522456427559e-5) * x 
					+ 3.03274662192286e-4) * e + w45 * ldata->weight[0];
				ldata->weight[2] = ((x * 1.34547929260279e-5 - 4.19389884772726e-4) * x 
					+ .00387706687610809) * e + w35 * ldata->weight[0];
				ldata->weight[1] = ((x * 2.09539509123135e-5 - 6.87646614786982e-4) * x 
					+ .00668743788585688) * e + w25 * ldata->weight[0];
				ldata->weight[0] = ldata->weight[0] - ldata->weight[1] - ldata->weight[2] -
					ldata->weight[3] - ldata->weight[4];
			} else {
				ldata->root[0] = r15 / (x - r15);
				ldata->root[1] = r25 / (x - r25);
				ldata->root[2] = r35 / (x - r35);
				ldata->root[3] = r45 / (x - r45);
				ldata->root[4] = r55 / (x - r55);
				ldata->weight[1] = w25 * ldata->weight[0];
				ldata->weight[2] = w35 * ldata->weight[0];
				ldata->weight[3] = w45 * ldata->weight[0];
				ldata->weight[4] = w55 * ldata->weight[0];
				ldata->weight[0] = ldata->weight[0] - ldata->weight[1] - ldata->weight[2] -
					ldata->weight[3] - ldata->weight[4];
			}
		}
	}
}

void Root6(Root * ldata) {
    double ff[19];
    double a[10], c[100]	/* was [10][10] */;
    long i, j, k, m, n,nn;
    double r[81]	/* was [9][9] */, s[100]	/* was [10][10] */, w[81]	/* was [9][9] */, x;
    long j1, k1, n1, jmax;
    double rt[10], dum;
    double root, poly, wsum;

/*	ITH ROOT OF THE JTH RYS POLYNOMIAL IS RETURNED IN R(I,J) WITH */
/*	THE CORRESPONDING WEIGHT FACTOR IN W(I,J).   J=1,2,...,N */
/*	THIS VERSION USES CHRISTOFFEL FORMULA FOR WEIGHTS, AND IT IS */
/*	COMPLETELY GENERAL UP TO NROOTS=9.  SEE THE REFERENCE */
/*		"NUMERICAL INTEGRATION USING RYS POLYNOMIALS" */
/*				H.F.KING AND M.DUPUIS */
/*			J.COMPUT.PHYS. 21, 144-165(1976) */

    n = ldata->nRoots;
    x = ldata->x;
    if (n < 2) n = 2;
    n1 = n + 1;
    nn = n + n;
    RysFun(&x, nn, ff);

    for (i = 1; i <= n1; ++i) {
		for (j = 1; j <= n1; ++j) {
			s[i + j * 10 - 11] = ff[i + j - 2];
		}
	}
    RysSmt(c, s, n1);

	for (i = 1; i <= n; ++i) {
		for (j = 1; j <= i; ++j) {
			w[i + j * 9 - 10] = 0.;
			r[i + j * 9 - 10] = 0.;
		}
	}

	wsum = ff[0];
	w[0] = wsum;
	r[0] = ff[1] / wsum;
	dum = sqrt(c[21] * c[21] - c[20] * 4. * c[22]);
	r[9] = (-c[21] - dum) * .5 / c[22];
	r[10] = (-c[21] + dum) * .5 / c[22];
	if (n == 2) goto L70;
	for (i = 3; i <= n1; ++i) {
		rt[i - 1] = 1.;
	}
	rt[0] = r[9];
	rt[1] = r[10];

	for (k = 3; k <= n; ++k) {
		k1 = k + 1;
		for (i = 1; i <= k1; ++i) {
			a[i - 1] = c[i + k1 * 10 - 11];
		}
		RysNod(a, rt, k);
		for (i = 1; i <= k; ++i) {
			r[i + k * 9 - 10] = rt[i - 1];
		}
    }

L70:
	for (k = 2; k <= n; ++k) {
		jmax = k - 1;
		for (i = 1; i <= k; ++i) {
			root = r[i + k * 9 - 10];
			dum = 1. / ff[0];
			for (j = 1; j <= jmax; ++j) {
				j1 = j + 1;
				poly = c[j1 + j1 * 10 - 11];
				for (m = 1; m <= j; ++m) {
					poly = poly * root + c[j1 - m + j1 * 10 - 11];
				}
				dum += poly * poly;
			}
			w[i + k * 9 - 10] = 1. / dum;
		}
	}

	for (k = 1; k <= n; ++k) {
		dum = r[k + n * 9 - 10];
		ldata->root[k - 1] = dum / (1. - dum);
		ldata->weight[k - 1] = w[k + n * 9 - 10];
	}
}

void RysFun(double * x, long n, double * ff) {
    double a, e;
    long i, m;
    double s, t, xx, fac, tol, sum, term, tmax, esave, facmin;


    /* Parameter adjustments */
    --ff;

	tol = 1e-14;
	xx = *x + *x;
	facmin = xx;

/*	VALUES X.GT.3200 AND NROOTS.EQ.6 GENERATE NEGATIVE ARGUMENTS */
/*	TO THE SQRT FUNCTION IN RYSSMT WHEN ONE USES E=10**(-35). */
/*	THE PROBLEM GOES AWAY, AND THE ROOTS AND WEIGHTS STABILIZE */
/*	FOR X OUT TO 5000 WHEN 10**(-40) TO 10**(-60) IS USED. */
/*	SINCE VAX UNDERFLOWS AT AROUND 10**(-38) WE HAVE TO BE MORE */
/*	CAREFUL THAN JUST SETTING E A FEW POWERS OF TEN DOWNWARD. */
/*	MICHEL'S ORIGINAL CODE IS SHOWN FIRST.   THE VALUE OF */
/*	360 IS WHERE IBMS GET INTO UNDERFLOWS (ABOUT 10**(-76)). */
/*	CHANGES INCLUDE THE COMPUTATION OF E IN THE ASYMPTOTIC */
/*	 CODE A FEW LINES DOWN. */

/* --   E     = 1.0D-35 */
/* --   IF(FACMIN.LT.360.0D+00) E=EXP(-X) */
/* --   IF(FACMIN.GT.80.0D+00) GO TO 100 */

	if (facmin < 360.) {
		e = exp(-(*x));
	} else {
		e = 0.;
    }
	if (facmin > 80.) {
		goto L100;
	}

	term = 1.;
	sum = 1.;
	fac = (double) (n);
	fac += .5;

L10:
	fac += 1.;
	term = term * *x / fac;
	sum += term;
	if (fac <= facmin) goto L10;
	t = term;
	s = sum;
	if (t > s * tol) goto L10;

	fac = (double) (n + n + 1);
	ff[n + 1] = sum * e / fac;
	m = n - 1;
	fac = (double) (m + m + 1);

L20:
	if (m < 0) return;
	ff[m + 1] = (e + xx * ff[m + 2]) / fac;
	--m;
	fac += -2.;
	goto L20;

/*     USE ASYMPTOTIC EXPANSION FOR LARGE ARGUMENTS. */

L100:
	if (e == 0.) {
		esave = 1e-35;
		for (i = 1; i <= 10; ++i) {
			e = esave * .1;
			if (e == 0.) goto L6;
		    esave = e;
		}
	L6:
		e = esave;
	}
	a = sqrt(.7853981633974483096156608 / *x);
	tmax = a * tol / e;
	term = 1. / xx;
	sum = term;
	fac = 1.;

L110:
	fac += -2.;
	term = fac * term / xx;
	sum = term + sum;
	t = term;
	if (fabs(t) > tmax) goto L110;

	ff[1] = a - e * sum;
	fac = -1.;
	m = 0;

L120:
	if (m == n) return;
	++m;
	fac += 2.;
	ff[m + 1] = (fac * ff[m] - e) / xx;
	goto L120;
}

void RysNod(double * a, double * rt, long k) {
    long k1, i, m;
    double p1, p2, r1, r2, p5, p6, r5, r6, r3, p3, r4, p4, dr, tol;
    double prod, delta, r;

/*	RETURNS IN RT(I) THE ITH ROOT OF A POLYNOMIAL OF ORDER K */
/*	WHOSE MTH COEFFICIENT IS STORED IN A(M+1).  IT IS ASSUMED */
/*	THAT THE INITIAL VALUES IN RT BRACKET THE FINAL VALUES. */

    /* Parameter adjustments */
    --rt;
    --a;

    /* Function Body */
    tol = 1e-11;
    k1 = k + 1;
    r2 = 0.;
    p2 = a[1];
    for (m = 1; m <= k; ++m) {
		r1 = r2;
		p1 = p2;
		r2 = rt[m];
		p2 = a[k1];
		for (i = 1; i <= k; ++i) {
			p2 = p2 * r2 + a[k1 - i];
		}
		prod = p1 * p2;
		if (prod >= 0.) {
			MessageAlert("RysNod: error.");
		}
		r5 = r1;
		p5 = p1;
		r6 = r2;
		p6 = p2;
	L30:
		r3 = r5;
		p3 = p5;
		r4 = r6;
		p4 = p6;
		r = (r3 * p4 - r4 * p3) / (p4 - p3);
		dr = r4 - r3;
		delta = dr;
		if (fabs(delta) < tol) {
		    goto L90;
		}
		dr *= .0625;
		r5 = r - dr;
		if (r5 < r3) {
		    r5 = r3;
		}
		r6 = r + dr;
		if (r6 > r4) {
		    r6 = r4;
		}
		p5 = a[k1];
		p6 = p5;
		for (i = 1; i <= k; ++i) {
			p5 = p5 * r5 + a[k1 - i];
			p6 = p6 * r6 + a[k1 - i];
		}
	L45:
		prod = p5 * p6;
		if (prod < 0.) {
		    goto L30;
		}
		prod = p3 * p5;
		if (prod > 0.) {
		    goto L60;
		}
		r5 = r3 * .25 + r5 * .75;
		p5 = a[k1];
		for (i = 1; i <= k; ++i) {
			p5 = p5 * r5 + a[k1 - i];
		}
		goto L45;
	L60:
		r6 = r4 * .25 + r6 * .75;
		p6 = a[k1];
		for (i = 1; i <= k; ++i) {
			p6 = p6 * r6 + a[k1 - i];
		}
		goto L45;
	L90:
		rt[m] = r;
    }
}

void RysSmt(double * c, double * s, long n) {
    long i, j, k, kmax;
    double v[10], y[10], fac, dot;

/*     ROUTINE RETURNS AN N BY N TRIANGULAR MATRIX C SUCH THAT */
/*     C(TRANSPOSE)SC=I,  WHERE I IS AN N BY N IDENTITY MATRIX. */

	for (i = 0; i < n; ++i) {
		for (j = 0; j <= i; ++j) {
			c[i + j * 10] = 0.;
		}
	}

	for (j = 0; j < n; ++j) {
		kmax = j;
		fac = s[j + j * 10];
		if (kmax != 0) {
			for (k = 0; k < kmax; ++k) {
				v[k] = 0.;
				y[k] = s[k + j * 10];
			}
			for (k = 0; k < kmax; ++k) {
				dot = 0.;
				for (i = 0; i <= k; ++i) {
					dot += c[i + k * 10] * y[i];
				}
				for (i = 0; i <= k; ++i) {
					v[i] -= dot * c[i + k * 10];
				}
				fac -= dot * dot;
			}
		}
/*		THE ORIGINAL CODE LEADS TO DIVIDE BY */
/*		ZERO ON VAX COMPUTERS FOR X>3100 OR SO. */
/* --		FAC=ONE/SQRT(FAC) */

		if (fac > 0.) {
			fac = 1. / sqrt(fac);
		} else {
			fac = 1e25;
		}

		c[j + j * 10] = fac;
		if (kmax == 0) continue;
		for (k = 0; k < kmax; ++k)
			c[k + j * 10] = fac * v[k];
    }
}
